#include <catch2/catch_test_macros.hpp>
#include "BVP.h"
#include <cmath>

using std::vector;

TEST_CASE("regular domian,u = x^2+y^2 , a = 1 , b = 0", "[dirichlet]")
{
    double epsilon = 1e-12;
    auto g = [](Point2d p)
    { return norm(p) * norm(p); };
    auto f = [](Point2d)
    { return -4.0; };
    BVP2d tmp({0, 0}, 0, f, BoundaryCondition(1, 0, g));
    vector solution = tmp.Solve(8);

    for (int i = 0; i < solution.size(); ++i)
    {
        REQUIRE(fabs(solution[i] - pow(norm(tmp.getGrid()[i]), 2)) < epsilon);
        // std::cout << tmp.getGrid()[i] << " : " << solution[i] << std::endl;
    }
}

TEST_CASE("regular domain,u = x^2+y^2 , a = 0 , b = 1", "[Neumann]")
{
    double epsilon = 1e-12;
    auto g = [](Point2d p)
    {
        if (max_of(p) == 1)
            return -2.0;
        else
            return 0.0;
    };
    auto f = [](Point2d)
    { return -4.0; };
    BVP2d tmp({0, 0}, 0, f, BoundaryCondition(0, 1, g));
    vector solution = tmp.Solve(8);

    double c = solution[0] - pow(norm(tmp.getGrid()[0]), 2);
    for (int i = 0; i < solution.size(); ++i)
    {
        REQUIRE(fabs(solution[i] - pow(norm(tmp.getGrid()[i]), 2) - c) < epsilon);
        // std::cout << tmp.getGrid()[i] << " : " << solution[i] - c << std::endl;
    }
}

TEST_CASE("regular domain,u = x^2+y^2 , a = 1 , b = 1", "[Mixed]")
{
    double epsilon = 1e-12;
    auto g = [](Point2d p)
    {
        if (max_of(p) == 1)
            return -2.0 + pow(norm(p), 2);
        else
            return 0.0 + pow(norm(p), 2);
    };
    auto f = [](Point2d)
    { return -4.0; };
    BVP2d tmp({0, 0}, 0, f, BoundaryCondition(1, 1, g));
    vector solution = tmp.Solve(8);

    // double c = solution[0] - pow(norm(tmp.getGrid()[0]), 2);
    for (int i = 0; i < solution.size(); ++i)
    {
        REQUIRE(fabs(solution[i] - pow(norm(tmp.getGrid()[i]), 2)) < epsilon);
        // std::cout << tmp.getGrid()[i] << " : " << solution[i] << std::endl;
    }
}

TEST_CASE("remove a disk centered at origin r = 0.5,a=1,b=0", "[Dirichlet]")
{
    double epsilon = 1e-12;
    auto g = [](Point2d p)
    { return norm(p) * norm(p); };
    auto f = [](Point2d)
    { return -4.0; };
    BVP2d tmp({0, 0}, 0.4, f, BoundaryCondition(1, 0, g));
    vector solution = tmp.Solve(8);

    for (int i = 0; i < solution.size(); ++i)
    {
        REQUIRE(fabs(solution[i] - pow(norm(tmp.getGrid()[i]), 2)) < epsilon);
        // std::cout << tmp.getGrid()[i] << " : " << solution[i] << std::endl;
    }
}

TEST_CASE("remove a disk centered at origin r = 0.4,a=0,b=1", "[Neumann]")
{
    double epsilon = 1e-12;
    auto g = [](Point2d p)
    {
        if (max_of(p) == 1)
            return -2.0;
        else if (min_of(p) == 0)
            return 0.0;
        else
            return 0.8;
    };
    auto f = [](Point2d)
    { return -4.0; };
    BVP2d tmp({0, 0}, 0.4, f, BoundaryCondition(0, 1, g));

    std::cout << "Begin test for irregular Neumann conditions...." << std::endl;

    int grid_size = 8;
    while (grid_size <= 64)
    {
        vector<double> solution = tmp.Solve(grid_size);

        vector<double> u, r;

        double norm1 = 0, norm2 = 0, norm_inf = 0;

        double c = solution[0] - pow(norm(tmp.getGrid()[0]), 2);

        for (int i = 0; i < solution.size(); ++i)
        {
            u.push_back(pow(norm(tmp.getGrid()[i]), 2));
            r.push_back(solution[i] - u[i] - c);

            norm1 += fabs(r[i]);
            norm2 += r[i] * r[i];
            norm_inf = std::max(norm_inf, fabs(r[i]));
        }

        std::cout << "grid_size : " << grid_size << std::endl;
        std::cout << "1-norm : " << norm1 / grid_size << std::endl;
        std::cout << "2-norm : " << sqrt(norm2 / grid_size) << std::endl;
        std::cout << "infinity-norm : " << norm_inf << std::endl;
        std::cout << std::endl;

        grid_size *= 2;
    }
}

TEST_CASE("remove a disk centered at origin r = 0.4,a=1,b=1", "[Mixed]")
{
    double epsilon = 1e-12;
    auto g = [](Point2d p)
    {
        if (max_of(p) == 1)
            return -2.0 + pow(norm(p), 2);
        else if (min_of(p) == 0)
            return 0.0 + pow(norm(p), 2);
        else
            return 0.8 + pow(norm(p), 2);
    };
    auto f = [](Point2d)
    { return -4.0; };
    BVP2d tmp({0, 0}, 0.4, f, BoundaryCondition(1, 1, g));

    std::cout << "Begin test for irregular Mixed conditions...." << std::endl;

    int grid_size = 8;
    while (grid_size <= 64)
    {
        vector<double> solution = tmp.Solve(grid_size);

        vector<double> u, r;

        double norm1 = 0, norm2 = 0, norm_inf = 0;

        for (int i = 0; i < solution.size(); ++i)
        {
            u.push_back(pow(norm(tmp.getGrid()[i]), 2));
            r.push_back(solution[i] - u[i]);

            norm1 += fabs(r[i]);
            norm2 += r[i] * r[i];
            norm_inf = std::max(norm_inf, fabs(r[i]));
        }

        std::cout << "grid_size : " << grid_size << std::endl;
        std::cout << "1-norm : " << norm1 / grid_size << std::endl;
        std::cout << "2-norm : " << sqrt(norm2 / grid_size) << std::endl;
        std::cout << "infinity-norm : " << norm_inf << std::endl;
        std::cout << std::endl;

        grid_size *= 2;
    }
}
